In [88]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import signal

When Models hurt prediction : a DC motor trial (p.132/155)

Définition du modèle et simulation

Modèle sous forme d'espace d'état :

$\begin{cases} \frac{d}{dt}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix} & =\begin{bmatrix}-\frac{R}{L} & -\frac{K}{L}\\ \frac{K}{J} & -\frac{f}{J} \end{bmatrix}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix}+\begin{bmatrix}\frac{1}{L} & 0\\ 0 & -\frac{1}{J} \end{bmatrix}\begin{bmatrix}V_{S}(t)\\ \tau_{L}(t) \end{bmatrix}\\ \begin{bmatrix}i_{obs}(t)\\ \omega_{Robs}(t)\\ V_{Lobs}(t) \end{bmatrix} & =\begin{bmatrix}1 & 0\\ 0 & 1\\ -R & -K \end{bmatrix}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix}+\begin{bmatrix}0 & 0\\ 0 & 0\\ 1 & 0 \end{bmatrix}\begin{bmatrix}V_{S}(t)\\ \tau_{L}(t) \end{bmatrix} \end{cases}$


In [89]:
def heaviside(x):
    if x == 0:
        return 0.5

    return 0 if x < 0 else 1

Simulation sorties moteur


In [110]:
@interact
def _(x0 = slider(0,10,10/100.0),
    x1 = slider(0,1,1/100.0),
    sim_duration = slider(10,30,default=10,label="Simulation duration"),
    sim_nb_samples = slider(10,100, default=100,label="Samples number"),
    steptime = slider(0,10,default=5,label="Step time"),
    stepgain = slider(-5,5,default=2,label="Step gain"),
    R = slider(0,2,2/100.0,default=1),
    K = slider(0,0.02,0.02/100.0,default=0.01),
    L = slider(0.5/100,1,1/100.0,default=0.5),
    J = slider(0.01/100.0,0.02,0.02/100.0,default=0.01),
    f = slider(0,0.2,0.2/100.0,default=0.1),
    action=selector(["MCC outputs", "RRA\'s"],label="Simulation choice")):

    ##### Modèle oracle
    A = matrix([[ -R/L, -K/L ],
            [ K/J , -f/J ]])
    
    B = matrix([[ 1/L , 0    ],
                [ 0   , -1/J ]])

    C = matrix([[ 1 , 0  ],
                [ 0 , 1  ],
                [ -R, -K ]])

    D = matrix([[ 0, 0 ],
                [ 0, 0 ],
                [ 1, 0 ]])

    ##### Modèle défectueux
    Af = matrix([[ -R*1.3/L, -K/L+1 ],
            [ K/J , -f/J ]])
    
    Bf = matrix([[ 1.2/L , 0    ],
                [ 0   , -1/J ]])

    Cf = matrix([[ 1 , 0  ],
                [ 0 , 1  ],
                [ -R, -K ]])

    Df = matrix([[ 0, 0 ],
                [ 0, 0 ],
                [ 1, 0 ]])

    ss_mcc = signal.lti(A,B,C,D)
    ss_faulty_mcc = signal.lti(Af,Bf,Cf,Df)
    
    t = np.linspace(0,sim_duration,sim_nb_samples)
    u1 = [stepgain*heaviside(x-steptime) for x in np.nditer(t)]
    u2 = np.array(np.zeros_like(t))
    U=np.vstack((u1,u2)).transpose()
    X0 = (x0,x1)
    tout, Y, X = signal.lsim(ss_mcc, U, t, X0)
    tout, Yobs, Xobs = signal.lsim(ss_faulty_mcc, U, t, X0)
    
    # Nice examples of interactive calculus : https://wiki.sagemath.org/interact/calculus
    # Plot help  : http://sage.brandoncurtis.com/plotting.html
    # Plot help2 : http://doc.sagemath.org/html/en/reference/plotting/sage/plot/plot.html
    # Matplotlib examples : http://matplotlib.org/examples/index.html
    
    if action == "MCC outputs":
        f, axarr = plt.subplots(3, sharex=True, figsize=(10,6))

        axarr[0].plot(t, Y[:,0], 'b')
        axarr[0].grid(True)
        axarr[0].set_ylabel('$i_{obs}(t)$',fontsize=20)

        axarr[1].plot(t, Y[:,1], 'r')
        axarr[1].grid(True)
        axarr[1].set_ylabel('$\omega_{R_{obs}}(t)$',fontsize=20)

        axarr[2].plot(t, Y[:,2], 'g')
        axarr[2].grid(True)
        axarr[2].set_ylabel('$V_{L_{obs}}(t)$',fontsize=20)

        plt.xlabel('Time (s)',fontsize=16)

        f.subplots_adjust(hspace=0.3)

        plt.show()
    elif action == "RRA\'s":
        n = len(t)-2
        U1 = np.zeros((n, 2))
        U2 = np.zeros((n, 2))
        Yobs1 = np.zeros((n, 3))
        Yobs2 = np.zeros((n, 3))
        Yobs3 = np.zeros((n, 3))
        r_ar1 = np.zeros(n)
        r_ar2 = np.zeros(n)
        r_ar3 = np.zeros(n)
        r_ir  = np.zeros(n)
        r_1 = np.zeros_like(t)
        r_2 = np.zeros_like(t)
        r_3 = np.zeros_like(t)
    
        ##### Calcul des RRA's auto-redondantes et inter-redondante
        for i in range(len(t)-2):
            U1[i] = U[i,:]
            U2[i] = U[i+1,:]
            Yobs1[i] = Y[i,:]
            Yobs2[i] = Y[i+1,:]
            Yobs3[i] = Y[i+2,:]

            ##### Calcul des valeurs des RRA's
            r_ar1[i] = J*L*Yobs3[i,0] - J*U2[i,0] - U1[i,0]*f + (K^2 + R*f)*Yobs1[i,0] + (J*R + L*f)*Yobs2[i,0] - K*U1[i,1]
            r_ar2[i] = J*L*Yobs3[i,1] - K*U1[i,0] + (K^2 + R*f)*Yobs1[i,1] + (J*R + L*f)*Yobs2[i,1] + R*U1[i,1] + L*U2[i,1]
            r_ar3[i] = K^2*Yobs1[i,2] + J*R*Yobs2[i,2] + J*L*Yobs3[i,2] + K^2*U1[i,0] + J*R*U2[i,0] - K*L*U2[i,1] + (R*Yobs1[i,2] + L*Yobs2[i,2] + R*U1[i,0])*f
            r_ir[i]  = J*L*Yobs2[i,2] + J*R*U1[i,0] - K*L*U1[i,1] + (J*R + L*f)*Yobs1[i,1] + (K^2*L + L*R*f)*Yobs1[i,0]

        ##### Calcul des RRA's directes
        for i in range(len(t)):
            r_1[i] = R*Yobs[i,0] + K*Yobs[i,1] + Yobs[i,2] - U[i,0]
            r_2[i] = Yobs[i,0] - Y[i,0]
            r_3[i] = Yobs[i,1] - Y[i,1]

        f, axarr = plt.subplots(7, sharex=True, figsize=(10,12))

        axarr[0].plot(t[0:-2], r_ar1, 'b')
        axarr[0].grid(True)
        axarr[0].set_ylabel('$r_{ar1}(t)$',fontsize=20)

        axarr[1].plot(t[0:-2], r_ar2, 'r')
        axarr[1].grid(True)
        axarr[1].set_ylabel('$r_{ar2}(t)$',fontsize=20)

        axarr[2].plot(t[0:-2], r_ar3, 'g')
        axarr[2].grid(True)
        axarr[2].set_ylabel('$r_{ar3}(t)$',fontsize=20)

        axarr[3].plot(t[0:-2], r_ir, 'g')
        axarr[3].grid(True)
        axarr[3].set_ylabel('$r_{ir}(t)$',fontsize=20)

        axarr[4].plot(t, r_1, 'b')
        axarr[4].grid(True)
        axarr[4].set_ylabel('$r_{1}(t)$',fontsize=20)

        axarr[5].plot(t, r_2, 'r')
        axarr[5].grid(True)
        axarr[5].set_ylabel('$r_{2}(t)$',fontsize=20)

        axarr[6].plot(t, r_3, 'g')
        axarr[6].grid(True)
        axarr[6].set_ylabel('$r_{3}(t)$',fontsize=20)

        plt.xlabel('Time (s)',fontsize=16)

        f.subplots_adjust(hspace=0.3)

        plt.show()



In [ ]: